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CALCULATION OF THE AERODYNAMIC LOADING OF SWEPT AND UNSWEPT FLEXIBLE 

WINGS OF ARBITRARY STIFFNESS 1 


By Franklin W. Diederick 


SUMMARY 


SYMBOLS 


A method is presented for calculating the aerodynamic load- 
ing, the divergence speed, and certain stability derivatives of 
swept and unswept wings and tanl surfaces of arbitrary stiffness. 
Provision is made for using either stiffness curves and root- 
rotation constants or structural influence coefficients in the 
analysis. Computing forms, tables of numerical constants 
required in the analysis, and an illustrative example are included 
to facilitate calculations by means of the method. 

INTRODUCTION 

The distribution of the aerodynamic loading on wings and 
tail surfaces is important both for the structural analysis of 
these components, since it determines the applied bending 
moment and torque acting at any station, and for their 
aerodynamic analysis, since it affects the stability deriva- 
tives to a large extent. At high speeds the aerodynamic 
loading, particularly in the case of swept wings, is greatly 
affected by the structural deformations caused by the load- 
ing. The present report is concerned with the determination 
of the effects of structural flexibility on the aerodynamic 
loading of wings of arbitrary plan form and stiffness. 

The present report treats the problem of aerodyn ami c 
loading by matrix methods. Aerodynamic induction is 
taken into account by means of approximate aerodynamic 
influence coefficients. When more accurate coefficients 
become available, they can readily be incorporated in this 
method. Structural flexibility is taken into account in the 
form of either calculated stiffness variations or measured 
influence coefficients. The required integrating matrices 
are presented for both a six-point and a ten-point solution. 
For the six-point solution convenient computing forms are 
included as well. The method is illustrated by means of 
an example. In addition to the analysis of the aerodynamic 
loading, the determination of the related divergence speed 
and of certain stability derivatives is discussed. 

For the convenience of the reader unfamiliar with matrix 
terminology, a summary of matrix methods has been includ- 
ed in the appendix. The sections entitled “Application of 
the Method” and, in particular, “Instructions for Solution” 
may be read without reference to the section entitled 
“Derivation of the Method.” 


A 


aspect ratio 



[A], [A*] aeroelastic matrices defined by equations (21) 
and (35) 

a section aerodynamic center, measured from leading 

edge, fraction of chord 
b wing span, inches 

b' wing span less fuselage width, inches 

c chord measured parallel to the air stream, inches 


c 


average wing chord, inches 




section lift coefficient ( — ) 

\qcf 

section lift-curve slope, per radian 

lift coefficient 

wing lift-curve slope, per radian 
effective lift-curve slope for twist distributions, 
per radian 



o 

[<?ll 


[C 2 ] 


[cy 


[cy 

EI 


GJ 

UJ 


root bending-moment coefficient 



T1 . . a- ■ . /Rolling moment\ 

rolling-moment coefficient qSb ) 

matrix relating concentrated and accumulated 
torques 

matrix relating concentrated loads and accumulated 
bending moments 

matrix converting torques due to distributed loads 
to torques due to concentrated torques 
matrix converting bending moments due to dis- 
tributed loads to bending moments due to con- 
centrated loads 

bending stiffness in planes perpendicular to the 
elastic axis, pound-inches 2 
location of elastic axis measured from leading edge, 
fraction of chord 

dimensionless distance along chord from reference 
axis to section aerodynamic center (e—a) 
torsional stiffness in planes perpendicular to the 
elastic axis, pound-inches 2 
integrating matrices for single integration from tip 
to root 

first row of matrix [I] 


1 Supersedes N"ACA TN 1876— Calculation of the Aerodynamic Loading of Flexible Wings of Arbitrary Plan Form and Stiffness, by Franklin W. Diederich, April 1949. 
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[IT], [IT] integrating matrices for double integration from 
tip to root 

[J/iJ first row of matrix [II] 

[I]" integrating matrix for single integration from root 
to tip (double transpose oLmatrix [I]) 

k dimensionless parameter ( -^} r - — tan a'') 

\(EI) r e x c r cos 2 A / 

wing lift-curve-slope ratio (C La jC L J 

lift, pounds 

running air load per unit length perpendicular to 
the plane of symmetry, pounds per inch 
accumulated bending moment (about axes parallel 
to the plane of symmetry unless specified other- 
wise), inch-pounds 
free-stream Mach number 
concentrated load, pounds 


K 

L 

l 

M 


Mo 

P 

pb 

2 V 

[Ql 

Qa i 0,9 

Qv 

g 


2 

8 

s 

T 


w 

Wc 


y 

a 

a 

r 

r 


) 


wing-tip helix angle 

aerodynamic influence-coefiicient matrix 
root-twist constants (see equation (15)) 
root-bending constant (see equation (15)) 
dynamic pressure, pounds per square inch 
dimensionless dynamic pressure 
( OL a q(b'/2ye lr c r i cos 
\ {GJ)r / 

reduced dynamic pressure (c^q tt 

dimensionless dynamic pressure 

( C La q(b'/2yc r tan A\ 

\ (j EI)r COS A / 

total wing area including part of wing covered by 
fuselage, square inches 

distance from wing root along reference axis, inches 
accumulated torque (about axes perpendicular to 
the plane of symmetry unless specified otherwise), 
inch-pounds 

concentrated torque, inch-pounds 
running torque due to air load about axes per- 
pendicular to the plane of symmetry, inch-pounds 
per inch 

fuselage width, inches 

distance between the effective root and the inner- 
most complete section of the torsion box per- 
pendicular to the elastic axis, inches 
lateral ordinate measured from plane of symmetry, 
inches 

lateral center of pressure, inches 
angle of attack, radians 

/ 6 \ 

equivalent angle of attack, radians + -q - 1 1 

local dihedral angle due to deformation, or slope of 
wing deflection curve at reference axis, radians 
structural deflection, inches 
lateral distance from wing root, inches 


A angle oLsweepback (measured to the reference axis 

unless specified otherwise), degrees 
influence-coefficient matrix for wing twist in pianos 
parallel to the air stream due to concentrated 
unit loads applied at the reference axis, radians 
per pound 

[<f> r ] influence-coefficient matrix for wing twist in planes 
parallel to the air stream due to concentrated 
unit torques applied in planes parallel to the air 
stream, radians per inch-pound 
cp angle of twist in planes perpendicular to the 

reference axis, radians 

Subscripts: 

c/2 midchord 

D _ .divergence 

jw flexible wing 

g geometric 

LE ' leading edge 

M due to bending moment 

MAC . pertaining to mean aerodynamic chord 

P due to concentrated load 

p damping in roll 

r at root or effective root unless specified otherwise 

rw rigid wing 

s structural (due to structural deformations) 

sub subsonic 

spr supersonic 

T due to torque 

TE trailing edge 

w wing exclusive of fuselage 

A in or pertaining to sections perpendicular to the 

reference axis 

Matrix notation: 

{ } . column matrix 

[ j row matrix 

[ j square matrix 

[ ] diagonal matrix 

[ ]' transpose of a matrix 

[ ]" double transpose of a matrix: first about tho princi- 

pal diagonal, then about the other diagonal 
[1] unit matrix 

[1/] matrix defined by equation (18) 

DERIVATION OF THE METHOD 

METHOD EMPLOYING STIFFNESS CURVES 

Assumptions. — In the development of the method the 
following assumptions are made: 

(a) All deflections and angles of attack are small. 

(b) The wing is mounted flexibly at an effective root per- 
pendicular to the elastic axis tlirough the intersection of the 
elastic axis and the fuselage (see fig. 1), the root rotations 
being proportional to the root bending moment and root 
torque. 
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Fcvubje 1.— Definition of geometrical parameters used in the analysis. 


(c) Ah elastic axis exists in the outer portion of the wing, 
this axis being defined as the elastic axis the wing would have 
if it were mounted rigidly some distance outboard of the root 
approximately perpendicular to the midchord line. (Near 
the root the elastic axis is defined as the extension of the out- 
board elastic axis.) 

(d) All deformations other than those due to the root 
rotations are given by the elementary theories of bending 
and of torsion about the reference axis, which in this case is 
the elastic axis. 

Air loads. — The force on a wing section of unit width 
parallel to the direction of flight is 

l=g_cc t (1) 

or, in matrix notation and in terms of the loading coefficient 
CCi 
Cr’ 

{Z}= 2 c r {^| (2) 


CCi 

The loading coefficient — - at any point on the span can be 

c T 

expressed in terms of the angles of attack at various points 
on the span by means of aerodynamic influence coefficients. 
For subsonic speeds, approximate aerodynamic influence 
coefficients may be calculated by the method of reference 1. 
With the resulting influence-coefficient matrix [Q\, the load- 
ing coefficients are given by the relation 

{^\=ClJQU«} (3) 


may be calculated more simply (but less accurately) by 
means of a modified strip theory; the section-lift distribution 
is rounded at the wing tips and reduced over the entire span 
by a factor which differs for angles of attack due to attitude 
and for angles due to any type of twist. On the basis of this 
approximation 

Cl = C LxCXg+C Lct a, "1 




(4a) 


=C. 


L a d 


for geometrical angles of attack which consist only of the 
angle of attack due to attitude, and 

Ci=C La (a, -fee.) ^ 

h (4b) 

J 


for all other geometrical angles of attack (due to rolling, 
or due to built-in twist, for instance). In equations (4a) 
and (4b), C La is an effective lift-curve slope for twist dis- 
tributions and a is an effective angle of attack defined by 


where 


a= or g + tea. 



Approximate values of C La and C La may be obtained by the 
reasoning of references 2 and 3 from the following equations: 


n A cos A 

A+2 cos A 

^ A cos A 

c La -c la CQS A 

so that 

A 4-2 cos A 
k ~A+ 4 cos A 


(5a) 

(5b) 

(5c) 


An influence-coefficient matrix which relates the loading 
coefficients and angle-of-attack values in the manner indi- 
cated in equation (3) may be obtained from equations (4a) 
and (4b): 

[Q]= |jf| [(i-‘)[u+4n] 

where [1J is a matrix defined by 


where a is the total angle of attack, which consists of the sum 
of the geometrical and structural angles of attack a t and a,. 
The geometrical angle of attack is that due to airplane alti- 
tude and built-in twist, whereas the structural angle of attack 
is that due to structural deformation. 

Instead of using the method of reference 1 the matrix [Q[ 


[ 1:1 = 


1 0 0 0 .. 

1 0 0 0 .. 

1 0 0 0 .. 

1 0 0 0 .. 
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CCi 

provided that in the column matrices of a and of — — . the 

values at the wing root or plane of symmetry are written at 

the top^. However, if the modified strip theory is used, 

equations (4a) and (4b) may be used directly without resort- 
ing to an influence-coefficient matrix so that, from equation 
(4a), 

{fH-Ul 131 (6a) 

and, from equation (4b), 

' ' (6b) 

For supersonic speeds the. loading coefficient given by 
unmodified strip theory can he expressed in the form 

- (6c) 

where c ta is the lift-curve slope pertinent to sections perpen- 
dicular to the midchord line. For instance, if the two- 
dimensional lift-curve slope is given by the Aclceret relation 

4_ 

then the section lift-curve slope to be used in equation (6c) is 

4 cos Ac /2 

*“ Yitfo 2 cos 2 Ae/ 2 — 1 

In the present report equation (6a) is used for the sake of 
definiteness. The modifications required in the method to 
use equations (6b) and (6c) rather than equation (6a) are 
obvious. 

If more accurate aerodynamic influence coefficients than 
those which correspond to the modified strip theory are 
desired, the coefficients of reference 1 .may be used in the 
manner described therein. Aerodynamic influence coeffi- 
cients for subsonic speeds can also be .obtained from the 
theoretical methods of references 4, 5, and 6 for calculating 
spanwise lift distributions. Whether the increase in accuracy 
obtainable by combining these methods with that of the 
present report warrants the corresponding increase in effort 

is somewhat questionable. _ . _ . . . 

Equations (2) and (6a) maybe combined to yield the follow- 
ing relation : 

{Z}=^ a 2C r [£]{«} (7) 

Similarly, the running torque t may be obtained from the 
relation 

t—G\d 

so that 

m=c^,vL^(£T|ra C8> 


(For cambered sections the pitching moment at zero lift must 
be added and the analysis of the following paragraphs modi- 
fied accordingly.) 

The accumulated torque T is obtained from the running 
torque and the running load by integrations inboard from the 
tip. The integration of the running torque may be performed 
by a matrix [/'] which is based on Simpson’s rule with a modi- 
fication suggested by V. M. Falkner at the lip. (See 
appendix.) The effect of Falkner’s modification is to round 
off the calculated load distribution and cause it to go to zero 
with an infinite slope at the tip, as the aerodynamic lift- 
distributions at subsonic speeds actually do. For supersonic 
speeds a similar matrix without lip modification [I] may be 
used. Both matrices are given in table I. The contribution 
of the running load to the accumulated torque is equal to the 
product of — tan A and the accumulated bending moment. 

The accumulated bending moment M is obtained by a 
double integration of the running load. The double integra- 
tion inboard from the tip may be performed by means of the 
matrices [77] and [II'] (given in table II), which are based on 
the equivalent of Simpson’s rule for moments, Palkner’s 
modification again being made at the tip in the case of [II']. 
The derivation of the integrating matrices is discussed in 
the appendix. 

In the following paragraplis subsonic flow will be assumed 
for the sake of definiteness. For supersonic flow, matrices 
[7] and [77] are used instead of [7'] and [II'\. The accumu- 
lated torque and bending moment may then be written as 

{T]=^[I'] {£}— tan A{il/J (9) 

and 

{M}=(0[77']{f} (10) 

The bending moments and torques referred to the elastic 
or reference axis, Ms and Ts, can be obtained from those 
given in equations (9) and (10) by means of the relations 

{Ma} =eos \{M] — sin A{T} 

sin A cos A (tttj [/'] Lq ('./]]!«! Cl) 

{2a} = cos A{T}+sin A{A/} 

=cosA Tj- [7'){f} 

=C La qc r \ y cos A[7'] {«} (12) 

Structural deformations. — The equations of equilibrium 
of a deformed wing referred to the elastic axis are 

aw 
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El (14) 

dy 

These equations must be integrated outboard from the root 
to obtain <p and F. The integrations may be performed by 
a matrix [/]" (see table III and appendix). Unlike the 
previously mentioned integrations this one is along the 
reference axis — that is, with s- as the independent variable 
rather than y or ij, which are the independent variables for 
integrations perpendicular to the plane of symmetry. 

b'12 

Consequently, the distance gog - must be used in conjune- 

tion with the dimensionless matrix [I\". 

To the deformations obtained in this manner the rotations 
<p T and r r due to the deformations of the root triangle must 
be added. The root rotations may be expressed by four 
dimensionless constants: 


Q*t — 


w t /(GJ) r 


(15b > 

ft ^ r T^ k r (;L5 C \ 

Qr T- Wt f(EI), (15c) 

ft (-1 r J \ 

Qr “~W'l(ET) r (15d) 

which may be combined into two other constants 

tan A<Jr 0 008 A (15e) 

i?, ») cos A <I6f) 

w t being defined as in figure 1. The deformations may then 
be written as 

( 16 ) 

and 

?^!^ ^Qr J ,[l/]{Z , i }] (17) 

where the matrix [1/] is defined by 

“o 0 0 0 . 

1 0 0 0. . 

1 0 0 0. . 

[1/1= 1 0 0 0 . . (18) 

1 0 0 0. . 


The angle of attack due to the structural deformations 
a, is related to <? and r by 

a,=(<p — r tan A) cos A (19) 

The aeroelastic equation. — If equations (11), (12), (16), 
and (17) are substituted in the matrix equivalent of equa- 
tion (19), the following relation is obtained: 

{ a ,}=g*Ul{«} (20) 

where the aeroelastic matrix [/lj is defined by 

[A] =[[/]" \j$]+yj 2 (G«r- Q-X tan A)[l/]+ 

WT, 1/1 U©T 

[*w L ttYw m [jfl 

( 21 ) 

the dimensionless dynamic pressure q*, by , 

* ^,ff(672) t «i r er I cos A 
2 ~ {GJ) t (22) 

and the parameter k, by 

k=@$- ^tt tan A (23) 

{EI) r ei T c r cos 2 A 

(The parameters g* and k are similar to the parameters a and 

— of reference 2.) 
a 

Solution of the aeroelastic equation. — If it is desired to 
calculate the aerodynamic loading corresponding to a given 
geometrical angle-of-attack distribution and dynamic pres- 
sure, equation (20) may be rewritten as follows: 

Ml- *g *M]J {«} = {*} (24) 

In this form it constitutes a set of linear simultaneous 
equations for values of a in terms of values of a t . From the 
calculated values of a the lift distribution may be deter- 
mined from equations (1) and (4). 

The divergence dynamic pressure may be obtained from 
equation (24) by setting the determinant of the square 
matrix on the left side of the equation equal to zero. This 
procedure is equivalent- to setting a g equal to zero in the 
term a of equation (20), so that 

{ a .} = Kq*[A}{*,} (25) 

The critical value of Kg* is then determined by matrix 
iteration and hence the divergence dynamic pressure from 
equations (5c) and (22). 

METHOD EMPLOYING INFLUENCE COEFFICIENTS 

The assumptions made in the preceding sections concern- 
ing the behavior of the wing structure are unnecessary if 
influence coefficients for the given structure are available 
from test data or refined methods of calculation. The co- 
efficients most convenient for this analysis are those giving 
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the rotation of the structure in planes parallel to the direction 
of flight due to vertical loads applied along a convenient 
reference axis and due to torques about lines perpendicular 
to the direction of flight. 

Since it is usually more convenient to apply concentrated 
rather than distributed loads in structural tests, the influence 
coefficients are considered in this analysis to have been ob- 
tained in this manner. These coefficients must be con- 
verted to coefficients which pertain to distributed loadings. 
This conversion may be accomplished either by means of the 
weighting matrices of reference 7 or by means of the conver- 
sion matrices described in this section. 

The angle of structural deformation a , may be expressed 
in terms of the influence coefficients and 4> P as follows: 

{ a ,}=[*r){T e }+[*rl{P} (26) 

where the T c ’s and P’s are arbitrary concentrated torques 
and loads, the latter being applied at the reference axis. The 
accumulated torques and bending moments about lines 
perpendicular and parallel, respectively, to the direction of 
flight may be related to the concentrated torques and loads 
by means of the summation matrices [Pi] and [Pa] (see 
appendix) as follows : 

{r} = [PJ{r < }— tan A{M} (27) 

{M}=^[P 2 ]{P} (28)* 

These relations may be solved for the values of T c and P 
required to produce given distributions of accumulated torque 
and bending moment 

{.TJ = [PJ-i({r}+tanA{M}J (29) 

{P}=y^[C 2 )-'{M} (30) 

The accumulated torques and bending moments produced 
by the air load are then 

[T}—~ [/'] {leie} — {M} tan A (31) 

{M}=(0[/P] {7} (32) 

Upon substituting equations (29), (30), (31), (32), (1), and (4) 
into equation (26), the following equation is obtained: 

= [^1 {«} (33) 

where 

2 (34) 

and 

[^=[e lr c r [<M [P 8 'l ICfi [£]] (35) 


where, in turn 

[p a ']=[p 1 r i [/ , i (36) 

[a , ] = [Pa]“ 1 [/P] (37) 

are given in tables IV and V. (Conversion matrices for 
supersonic flow [P 3 ] and [P 4 ] can be calculated from [/] and 
[II], if needed, by replacing [/'] and [//'] with [/] and [77] 
in equations (36) and (37).) 

The solution of equation (33) is obtained in the manner 
previously described for equation (20). 

APPLICATION OF THE METHOD 
DETERMINATION OF THE STRUCTURAL PARAMETERS 

At the time an aeroelastic analysis is performed no experi- 
mental stiffness data are usually available, so that either the 
calculated stiffness curves or calculated influence coefficients 
must be used. In order to use the stiffness curves it is neces- 
sary to assume the existence of a reasonably straight elastic 
axis. The location of this axis may be estimated by con- 
sidering it to be the line connecting the shear centers of the 
individual sections. If the elastic axis obtained in this 
manneris not reasonably straight within one or two percent 
of the chord, the results of the analysis may not be suffi- 
ciently reliable. . 

The stiffnesses GJ and El do not have much physical 
significance inboard of the last point where there is a com- 
plete cross section of the torsion box. (Seo fig. 1.) In order 
to arrive at estimates of the root stiffnesses (GJ) r and {EI)„ 
which serve primarily as reference values in this analysis, 
the stiffness curves have to be extended. It is convenient 
to consider the stiffnesses to be constant inboard of Lite last 
complete section of the torsion box; this procedure should 
yield conservative values of the root rotations. 

The most difficult problem incurred in analyzing the de- 
flections on the basis of stiffness curves appears to be the 
estimation of the root rotations. As used in this analysis, 
they are the torsion and bending deflections imposed by the 
triangular inner portion of the wing and the carry-through 
bay on the rest of the wing. As seen in figure 2, which is 
plotted from the data of reference 8, these values are essen- 
tially constant along the span, so that they actually consti- 
tute rigid-body rotations. (The bending rotations have 
been obtained by taking the difference in slopo between 
curves calculated by considering the wing to be cantilevered 
at the effective root — the root used to calculate torsional 
deformations in reference 8 — anti the averages of the leading- 
edge and trailing-edge deflections actually measured. The 
twists were obtained by subtracting the twists calculated on 
the basis of the assumed effective root from the measured 
twists.) 

The rotations should, in any practical case, be calculated 
by analyzing the triangular root and the carry-through bay 
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Rotation Av. value Q value 



Fiacre 2— Rotations of a 45° swept box beam due to root deflections (data from reference 8). 
ti=104 inches; w .= 15 inches; (GJ r ) r =5.WX10 a pound-inch'; (El) ,=9.47X10* pound-inch'; 
TV =43,420 inch-pounds; Af,=260,000 inch-pounds. 


The influence coefficients used in the analysis consist of 
the rotations of sections parallel to the direction, of flight 
due to concentrated unit torques in planes parallel to the 
plane of symmetry or concentrated unit loads at the refer- 
ence line. When measured, these rotations (in radians) 
may be entered in tables of the form: 


[*ri 

TWIST AT STATION ^ DIJE T0 UNIT CONCENTRATED 
TORQUE AT ^ 


i 

472 

0-2 

0.4 

0.6 

0.8 

0.0 

1-0 

0 







0.2 







0.4 







0.6 







0.8 







0.9 








and be made dimensionless by means of equations (15a) to 
(15f). If such an analysis is not available, the dimensionless 
rotation parameters shown in figure 2 may be used as a 
guide. It- must be kept in mind, however, that in the case 
of a sweptforward wing the parameters Q rw and Qr T would 
have the opposite sign and that for antisymmetries! loadings 
the root rotations may be slightly different than for sym- 
metrical loadings (see reference 9), 

If higher-order structural effects are to be taken into 
account — such as shear lag, bending stresses due to torsion, 
local stresses due to cut-outs, discontinuity of the elastic 
axis, and so on — structural influence coefficients may be 
calculated by means of the method of reference 10 and used 
in the present analysis in the same manner as measured 
influence coefficients. 

Once the structure under investigation is built, fairly 
simple deflection tests, similar to those performed in refer- 
ences 8 and 9, may be used to check the root-rotation 
parameters by calculating the differences between the 
observed rotations and those calculated by simple beam 
theory considering the wing cantilevered at the effective 
root; at the same time the existence and estimated location 
of the elastic axis may be verified. If the experimental 
program is fairly extensive it is desirable to measure influence 
coefficients directly. These influence coefficients can then 
be used in conjunction with the alternate method described 
in the preceding section to obtain a quick check on the aero- 
elastic analysis based on calculated stiffnesses. 


[**] 

TWIST AT STATION ^ DUE TO UNIT CONCENTRATED 
LOAD AT ^ 


472 

0.2 

0.4 

0.6 

0.8 

0.9 

1.0 

0 







0.2 







0.4 







0.S 







0.8 







0.9 


- 






These particular tables would be used for a six-point anal- 
ysis; similar tables would be used for a ten-point analysis. 
In either case it is to be noted that the twists are measured 

at values of from 0 to 0.9, whereas the loads are applied 


at values of 


b'[2 


from 0.2 to 1.0. The tables obtained in this 


maimer constitute the desired influence-coefficient matrices. 

If the wing sections are found to twist nonuniformly, so 
that they become cambered in effect, the angles of twist a, 
to be entered in the influence-coefficient matrices have to be 
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defined in a different manner according to whether the aero- 
elastic analysis is performed for subsonic or supersonic speeds. 
At subsonic speeds the lift depends on the slope of the mean 
camber line at the three-quarter-chord point, so that the 
effective angle of attack is 


a s 


_n (fc/2 — firs) 

x S- 


(38) 


At supersonic speeds the lift depends primarily on the aver- 
age slope of the mean camber line, so that 


T LB — £ TE 



(39) 


DETERMINATION OF THE AERODYNAMIC PARAMETERS 

The selection of the aerodynamic parameters C Lat and e , 
for the calculation of the divergence speed has been discussed 
in referenced. For calculating the aerodynamic loading at 
a given flight condition the aerodynamic parameters are 
chosen for that flight condition. The use of the effective lift- 
curve slopes 0 La and C La ^ in. conjunction with modified strip 
theory is applicable only to subcritical Mach number’s. At 
higher speeds no simple span correction is available; neglect 
of the span correction tends to be conservative for calculation 
of the divergence speed and the aerodynamic loading, 
however. 

INSTRUCTIONS FOR SOLUTION 

Two sets of integrating matrices have been prepared, one 
for a six-point solution and one for a ten-point solution. 
The former should be adequate for all practical purposes; 
only where the stiffness curves are very irregular near the 
root does the ten-point solution have to be used. The points 
considered by the two sets of tables in the case of subsonic 

flow are 0.2, 0.4, 0.6, 0.8, and 0.9 for the shorter 

solution and 1 ^ 2 = 0 , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, and 
0.9 for the longer solution. For supersonic cases an addi- 

yi 

tional point at JY2 =1 '° * s use< ^ ^ 01 ’ solutions. The 

procedure to be followed for all solutions is identical; although 
computing forms are presented in this report only for the 
six-point solution, their extension to the ten-point, solution 
is obvious. 

Calculation of the matrices. — The first step in the aero- 
elastic analysis by means of the stiffness curves is the calcu- 
lation of the aeroelastic matrix [/4] from the physical and 
geometrical parameters of the wing. These parameters are 
conveniently tabulated in a form of the type shown in table 
VI (a). The computation is then carried out according to 
the instructions of table VI (b), each, step in the procedure 
being identified by the number in the upper left corner of 
each box. It must be kept in mind that many of the 
operations call for matrix multiplications where the order of 


the multiplicands is of importance. (A brief summary of 
matrix methods is presented in the appendix.) The aero- 
elastic matrix is obtained as the last step (step 13) of the 
computations in this form, which constitutes an evaluation 
of the matrix [At] given in equation (21). 

In the computing form of table VI(b) and in the illustrative 
example the matrices II'] and [II') are used for the supersonic 
case. This procedure results in a slight saving in effort and 
is justified to a certain extent because even in supersonic flow 
the reduction of the actual lift carried by a section compared 
with the strip-theory value is largest at the tip. In general, 
however, use of the matrices [I] and [II] would probably be 
preferable for supersonic speeds. 

A special case arises when e lr is zero. If e i is not zero 
along the remainder of the span, its value at some point 
other than the root may be used as a reference value. The 


, . e, /e\n 

matrix — ( — 1 
L«i r W 


and the multiplying factors of steps S 


and 9 as well as the definition of the parameter q* arc then 
based on the value of at this other reference station. If e, is 
zero along the entire span, step 1 and. steps 3 to S may be 
omittecLand steps 9 to 13 should be modified as follows: 


Step 9 


(EI) r w, Q«,t M 
{GJ) r b'/2 tan A 1,1 


Step 10 [®]— [®1 


Step 11 As is 


Step 12 Omit 


Step 13 [-A] ei -o=[©l[@l 


If structural influence coefficients of the proper type are 
available, the calculation of the aeroelastic matrix [A f ] is 
carried out directly by means of equation (35). 

Solution for divergence dynamic pressure.— In order to 
determine the value of the dynamic pressure corresponding 
to divergence, the aeroelastic matrix [A] or [/l f ] is iteratod 
(see appendix) as indicated by equation (25). Table VII (a) 
may be used for this purpose. The result is the critical 
value of Kq* or nq*. The divergence dynamic pressure is 
then calculated from equations (22) or (34). (It is to bo 
noted that this pressure will be in pounds per square inch.) 
Since the aeroelastic matrix is independent of tho Mach 
number, except insofar as c t varies with Mach number, the 
same critical value of nq* or Kq} may be used to calculate 
the divergence dynamic pressure for an entire range of 
Mach numbers. If the value of e t changes, however, ns it 
does between the subsonic and supersonic region, critical 
values of q* and q} have to be calculated for both values 

Of 

If the value of e x is zero along the entire span and the 
matrix [A] has been calculated according to the modified 
instructions, iteration of the matrix will give the value of 
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the parameter Kg at divergence. From the definition of g 


*2 = 



tan A 


( ET)r COS A 


(40) 


the divergence dynamic pressure may then be calculated. 

Solution for aerodynamic loading. — In order to calculate 
the aerodynamic loading corresponding to a given flight 
condition and geometric angle-of-attack distribution, the 
aeroelastic matrix [A] or [A*] is multiplied by the value of Kg* 
or Kg* calculated for the given flight condition and subtracted 
from the unit matrix [l]. (See equation (24).) The result 
may be entered in table VII (b). Again it must be noted 
that the value of the aeroelastic matrix varies with the flight 
condition if e, varies, so that the aeroelastic matrix corre- 
sponding to the proper value of e t must be selected. 

The elements of the resulting matrix are the coefficients 
of a set of simultaneous linear algebraic equations for the 
unknown values of the effective angle-of-attack distribution 
of the deformed wing {« } in terms of the known angle-of- 
attack values of the rigid wing {a,}. Table VII (b) is set 
up for the calculation of the additional loading, the damping- 
in-roll loading, and a third arbitrary loading; as many 
loadings as desired may, of course, be calculated by this 
method. The solution of the equations may be carried out 
in any convenient manner. The form of table VII (b) has 
been prepared for use in conjunction with Crout’s method 
of solving linear simultaneous equations (reference XI). 

The values of a may also be obtained by means of an 
iterative rather than a simultaneous-equation type of solu- 
tion. A first approximation to the structural twist may be 
calculated by premultiplying the column of the values of the 
geometrical angle of attack a t by the matrix [A] and 
multiplying the resulting column by g*. In so doing, the 
contribution of a x to the total equivalent angle a has been neg- 
lected. However, the approximate values of a, calculated 
in this maimer may be multiplied by the factor k and added 
to the values of the geometrical angle of attack to obtain 
appro xim ate values of a. The column of these values is then 
pr em ultiplied by q*[A] to obtain a better approximation to a,. 
This procedure may be repeated until the solution converges. 

Alternatively, the first approximation to the column of a, 
may be multiplied by k and premultiplied by q*[A] to obtain 
approximately the contribution of the term a, in a to the 
calculation of a, by means of equation (20). In turn, the 
effect of this contribution on a, may be calculated by pre- 
multiplying the contribution by Kg*[A], and so on. The final 
value of a, is then the sum of all these contributions. 

Bo th of these procedures are equivalent. They may be sys- 
tematized and shortened as follows: The matrix [A] is 
entered at the upper left of table VH (c). Two sets of 
values of {a,} are listed at the left of the table. For the 
rolling case, 


y b' n ,w 

a ‘~bl2~b 6'/2 ‘ & 


( 41 ) 


These columns are then premultiplied by the matrix [A]. 
The resulting columns are entered in the columns to the right 
of the columns of {a t } and again premultiplied by the matrix 
[A]. The results of this second matrix multiplication are 
entered in the next column and premultiplied again by the 
matrix [A]. This procedure is repeated until all the ele- 
ments of one column, say the rth, differ from the elements of 
the preceding column, the (r-l)th, only by a constant factor 
(within 1 percent or less). This constant factor is equal to 


t — jr — and is entered at the bottom of the table. The 
(*2 *)d 

dimensionless dynamic pressure at divergence is, conse- 
quently, obtained automatically in this iterative procedure 
and need not be calculated separately by means of table 

VH (a). 

The dynamic pressures of interest and the corresponding 
values of Kg* are entered in the matrix at the upper right of 
table VTI(c). Also entered are the corresponding values of 
(/eg*) 2 , (xg*) 3 , and so on, up to and including (k ; g*) r_1 . In 

( K q*v 

the next row, however, the values of , - r — / — r are entered 

1 — (2/2/>J 


instead of (/cg*) r . 

The values of a for the two sets of values of a t and for 
the various values of g of interest are then obtained by pre- 
multiplying the matrix which consists of the rows Kg*, (/eg*) 2 , 
(sg*) 3 , . . . (including the row of I’s) by the matrices of the 
columns {a t }, [A]{o/J, [A] 2 {a t } ; and so on. The resulting col- 
umns comprise the desired values of a for the various cases. 
These values are entered in the appropriate columns of 
table VII (c). 

For most conventional plan forms and structures the 
results of the iterative procedure described in the foregoing 
paragraphs converge in two to four cycles (r=2 to 4); how- 
ever, for a wing which is uncommonly flexible near the tip, 
more cycles may be required. In general, the simultaneous- 
equation type of solution is preferable when several geo- 
metrical angle-of-attack conditions are to be analyzed at one 
or two dynamic pressures. The iterative procedure, on tne 
other hand, is preferable when only one or two angle-of- 
attack conditions are to be analyzed for several dynamic 
pressures. When many stations along the wing span are 
taken into account, the iterative procedure is preferable in 
almost all cases. 

In the case where e x is zero along the span, the headings 
at the top of table VH (b) should be modified to read 


[Lll— Kg[Al, 1=0 ] 

where [A] ej=0 has been calculated according to the modified 
instructions and Kg has been obtained by iterating [A] rt=0 - 
(See also equation (40).) 

The values of {«} calculated for the additional load dis- 
tribution (a s = 1) constitute values of the ratio c tf JCi ra or 
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(cc ; ) fwl(cci) TW in view of the assumptions made concerning the air 
forces. The section loading of the flexible wing is obtained 
from the relation 


CCi = cC L a 

Ct 

(42a) 

or in dimensionless form 


*1$ 

II 

R 

ftl 

(42b) 


The wing lift and the bending moment at the wing root may- 
then be obtained by integrating the load distribution. These 
integrations may be performed: conveniently by multiplying 
the values of cci/c r by the first rows of the matrices [l'\ and 
[IF], respectively. Thus 

L a =2<iCi, B c r ^ Ui'J [_£]{«} (43) 

and 

M t =qC> m Cr [£\ {«} (44) 

These quantities may be expressed as dimensionless coefficients : 

/f 

and 

The lateral center of. pressure of the wing load y may be 
determined from the relation 


V 


w 

2 



( 45 ) 


The fore-and-aft location of the aerodynamic center of the 
wing load measured rearward of the leading edge of the 
mean aerodynamic chord as a fraction of the mean aero- 
dynamic chord may then be estimated from the relation 


a | y y mac 
Cmac CmaC 


tan A a 


(46) 


where A a is the sweep of the section aerodynamic center line. 

For any other geometrical angle-of-attack distributions, 
such as those due to built-in twist or those due to rolling, 
the same section lift-curve slope should be used for the 
geometrical as for the structural deformations, so that Cl u 
is replaced by Ck„ e , and k is unity, and a in equations (42) 

to (44) is replaced by a. 

For the damping-in-roll distribution with a tip helix angle 
of 1 radian, a e is given by equation (41) and the rolling- 
moment coefficient due to the whig load is given by 


Gi ,—2 


M r 

qSb 


~ (47) 


where M r is obtained from equation (44) with the values of 
a calculated for this case. 


The contribution of the wing to other stability derivatives 
may be obtained, similarly, by integrating the load distribu- 
tions due to the angle-of-attack distributions caused by the 
motion of interest, as described in reference 12; in the case 
of swept wings, particular care must be taken in selecting the 
proper angle-of-attack distribution and in accounting for 
the lateral inclination of the lift vector. (See reference 3.) 

If the aerodynamic loading or the stability derivatives 
are to be obtained for a wide variety of flight conditions, it is 
convenient to systematize the calculations in the following 
manner.;. The aeroelastic matrix is computed for both the 
subsonic and supersonic aerodynamic-center values and 
iterated for both cases to obtain the subsonic and supersonic 
values .of the divergence parameter From these 

values the divergence dynamic pressure may be computed 
by means of equation (22) and plotted against Mach num- 
ber, as suggested in reference 2; on the same plot, values 
of the actual dynamic pressure may be plotted against 
Mach number for various altitudes of interest. Such a plot 
for a wing, the physical characteristics of which are given in 
figure 3, is shown in figure 4. 

Since at a given Mach number the ratio Kq*/(Kq*) D is 
equal to the ratio qJq D , the range of values of k2*/( k 2*)o of 
interest may be established from this plot for both the sub- 
sonic and the supersonic regions. Several representative 
values of this ratio may then be chosen within the given 
ranges and the corresponding values of kq* computed from 
the previously calculated values of («£*) D . The aerodynamic 
loading is calculated for these values of Kg* by using the 
appropriate matrix [ML] and may be plotted in the form of 
(cci)/J(cci) ra , with the ratio qjq D as a parameter. From 
these curves (or from the values of a) the wing lift coefficients 
may be obtained and plotted in the form (CO/J^Cz,)™ 
against g/gej the other coefficients may be obtained and 
plotted in a similar form. 

For any specific flight condition the value of qjq n may 
then be obtained from the plot of q and q D against Mach 
number. The loading, lift coefficient, or other item of 
interest-may be obtained from the plots which give these 
items in terms of the rigid-wing values. Once the rigid- 
wing values at the given Mach number are known, the 
flexible-wing values may then be obtained immediately". 

ILLUSTRATIVE EXAMPLE 

In order to illustrate the method described in the preced- 
ing sections, a swept wing has been analyzed. The physical 
and geometrical parameters of the wing are shown in figure 
3 and the upper part of table VIII (which follows the form 
of table VI (a)). The cliord, Ci c*, and stiffness matrices 
have been obtained from the given parameters and arc shown 
in the lower part of table VIII. 

The calculation of the aeroelastic matrix for the subsonic 
case has been carried out by means of the form of table \T(b). 
All but three of the steps of the computation are shown in 
table IX, numbered in the same order as in table VI(b). 
Steps 1, 2, 6, 7, 11, and 12 constitute matrix multiplications, 
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Fioc'ee 3.— Parameters of the example wing. 

which are carried out in the order indicated; steps 5 and 13 
constitute matrix additions or subtractions; steps 3 and 4 
constitute multiplications of matrices by constants. 

The aeroelastic matrix is iterated in table X (a) (which 
follows the form of table VII fa)) to yield a value of 
( K q*) D = -2.208 

From this value and a value of this parameter com- 
puted in the same manner for supersonic speeds (using 
C' L ~C h =C[ ), the divergence dynamic pressure has been 

CC OCp ot 

calculated by means of equation (22) on the basis of esti- 
mated values of the effective lift-curve slope. The variation 
with Mach number of the divergence dynamic pressure, the 
actual dynamic . pressure at sea level, and the estimated 
effective lift-curve slope is shown in figure 4. 

For a value of —=—0.25, such as would be obtained ap- 
2 o 

proximately at a Mach number of 1.0, "the aerodynamic 
loading has been calculated for the additional-angle-of- 
attack case and the damping-in-roll case in table X(b), 
which follows the form of table VTI (b) . The values of a s for 
the damping-in-roll case have been calculated from equa- 
tion (41). The aerodynamic loadings, in addition to those 
calculated for other values of q/q D , have been plotted in 
figure 5 as ratios of the flexible-wing loadings to the rigid- 
wing loadings. The curves have been integrated to yield 
wing lift and rolling-moment coefficients as well as the 
aerodynamic center of the wing load, which are shown in 

table X(b) for the case of ^-=—0.25 and which are plotted 
against —-~'m figure 6. 

930046 — 51 39 



Figure 4>— Effect of Mach number on the divergence dynamic pressure and lift-curve slope 
of the example wing. 


The wing lift coefficient is defined in such a manner 
that if the fuselage lift is known and made dimensionless by 
dividing by q and S the resulting fuselage lift coefficient may 
be added directly to the wing lift coefficient. This definition 
and the fact that figure 5 (a) is plotted over the fraction of 


the wing-alone span explain the fact that the area under 


the curve of figure 5 (a) is not 1. The aerodynamic center 
as plotted in figure 6 constitutes the center of pressure of 
only the wing load. In order to obtain the airplane aero- 
dynamic center, the magnitude and center of pressure of 
the fuselage load would have to be known and taken into 
account. 


DISCUSSION 


Both the aerodynamic and the structural assumptions 
made in this analysis are somewhat more realistic than those 
made in reference 2. The device employed in this analysis 
of calculating the air forces for wing sections parallel to the 
direction of flight- and then transferring them to sections 
perpendicular to the elastic axis obviates the necessity of 
replacing the actual wing with one the root and tip of which 
are perpendicular to the elastic axis for the purpose of 
analysis. 
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Figure 5.— Load distribution of example wing. 


The inclusion of Falkner’s modification (see appendix) 
in the integrating matrices has the effect of rounding off the 
load distribution approximately in the manner observed at 
subsonic speeds. At supersonic speeds the load distributions 
do not go to zero in the same manner, but even at supersonic 
speeds there is some reduction of load at the tip, the total 
magnitude of which is not far from the reduction obtained 
by lalkner’s modification. However, the use of special 
matrices [/] and [II] is indicated if the amount of this 
reduction is known. ... ... 

The assumption, that induction effects may be approxi- 
mated by an over-all reduction of the strip theory loading 
(rounded off as previously described) at subcritical speeds and 
may be neglected at supersonic speeds, may be avoided by 
using aerodynamic influence-coefficient matrices instead 
of the effective lift-curve slopes. Available methods of 
calculating such influence coefficients from theoretical 
methods for calculating lift distributions for wings of arbi- 
trary plan form at subsonic and supersonic speeds are gener- 
ally either too inaccurate or too time consuming for practical 
purposes. The empirical method of reference 1 has the 
advantage of simplicity and fairly good accuracy compared 
with theoretical methods but is applicable only to sym- 
metrical lift distributions. . . 

Although the analysis of this report has been performed 
for wings consisting of uncambered_sections, the analysis 
is directly applicable as well to the determination of the 
additional loading of wings with cambered sections. The 
loading of such wings due to the section pitching moment 
at zero lift may be determined by modifying the analysis 
somewhat. 




Figure 6.— Lift coefficient, rolling-moment coefficient, and aerodynamic center of eximplo 

wing. 

The assumption of an effective root perpendicular to the* 
elastic axis made in reference 2 for the purposes of calculating 
the structural response is carried over in this analysis. 1 1 is 
modified, however, to the extent that the root is no longer 
considered to be rigid as in reference 2, but flexible, both in 
torsion and bending. It has been demonstrated in reference 
8 that the deflections of a swept beam may bo estimated on 
that assumption, provided the root-rotation parameters are 
known. By assuming the effective root at the intersection 
of the elastic axis with the side of the fuselage, the root 
bending due to bending moment and root twist due to 
torque are minimized. The bending due to twist and twist 
due to bending are the same regardless of the location of the 
effective root. 

The method of introducing the root rotations into the 
analysis by means of the matrix [1/] assures that the struc- 
tural twist in planes parallel to the direction of flight is zero 
at the fuselage. In so doing the assumption, is made that- 
the part of the wing structure within the fuselage docs not 
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deform under load. From figure 2 it is seen that the local 
values of the root rotation either tend to approach zero at the 
root or tend to cancel each other. (The root rotations shown 
in figure 2 do not contain the rotation due to deformation of 
the cany-through bay; for cases in which this rotation is 
believed to be sizeable it has to be taken into account in the 
coefficients Q < and Qa r ) If the root-rotation constants 
are known, the structural deformations can therefore be 
predicted quite accurately by the assumptions made. 

The manner in which the equations of equilibrium are 
solved by means of the integrating matrices accounts for the 
true chord and stiffness variations. It does not necessitate 
replacement of the actual wing by constant-chord segments 
with all the flexibility concentrated at the ends of the seg- 
ments, an approach which has been used extensively in the 
w'ork on aeroelastic problems of unswept wings. 

A further refinement which obviates the necessity for 
making any structural assumptions other than that of small 
deflections is the use of measured influence coefficients or 
coefficients calculated by refined structural theories in the 
aeroelastic analysis. Wherever such coefficients are avail- 
able it is, of course, of advantage to use them. 

No explicit account has been taken in the analysis of the 
effects of the inertia loadings on the structural deformations 
and hence the aerodynamic loading. On swept wings, in 
particular, their effects may be considerable. For the pur- 
poses of this analysis, however, the structural deformations 
due to inertia loading may be considered part of the geo- 
metrical angle of attack and the rigid-wing geometrical angle 
of attack may be modified accordingly. The deformations 
due to the inertia loading may, incidentally, be calculated 
conveniently by means of the matrices [7], [77], and [7]". 

Some of the general observations made in reference 2 con- 
cerning the divergence phenomenon are corroborated by the 
example. As expected of a wing with a considerable amount 
of sveepbaek, the divergence dynamic pressure is negative; 
consequently, the wing cannot diverge. The divergence 
dynamic pressure is useful as a reference value, however; 
the values of the load distribution and the stability param- 
eters divided either by the corresponding rigid-wing values 
or by the section lift-curve slope depend only on the ratio 
of the actual to the divergence dynamic pressure. 

The type of plot shown in figure 4 is therefore quite useful 
in the analysis of aeroelastic phenomena. As pointed out in 
reference 2, this chart may also be used to estimate the 
actual divergence dynamic pressure where there is a possi- 
bility that the wing may diverge. It appears that the crit- 
ical values will tend to occur at either extremity of the 
transonic region. 

As would be expected qualitatively, the effect of wing 
flexibility in the case of the example wing is to unload the 
wing tips owing to the fact that they bend up. The lift- 
carried by the wing is therefore Less than that carried by a 
rigid wing, the center of pressure being farther inboard and 
the aerodynamic center being farther forward. 

The difference between the supersonic and subsonic values 
of the loading, the lift and rolling-moment coefficients, and 
the aerodynamic center for a given value of q{q D is due to 


the difference in the values of e u 
Another item of possible interest is the fact that the vari- 
ations of the parameters (/eg) D and (nq*) D (the parameters 
d D and a D oi reference 2) for the example problem are approx- 
imately linear (see fig. 7), as would be expected from the 
results of the analysis of reference 2. The deviations from 
linearity are most pronounced near the points for (£q) D =b 
(that is, A=0). They are due to the effects of the root 
rotations, in particular, the bending due to torsion^ and 
torsion due to bending; these effects were neglected in the 
approximate analysis of reference 2. The points of figure 7 
correspond to the example wing and the wings that would 
be obtained by rotating the example wing to the unswept 
and 37.5° sweptforward positions in such a manner as to 


e,c T cos 2 A (El) 

keep constant the parameters ' ( GJ) ’ aS as 

the chord, stiffness, and moment-arm distributions g,. Points 
are shown for both the subsonic and supersonic variations 
as well as for the case when g,=0 over the entire span 
(( K 2*)z>=0). The difference between the subsonic and super- 
sonic lines is due entirely to the difference in the distributions 
of g, rather than the difference in the magnitudes of e lr or 


in the character of the lift distributions. 

The present analysis is concerned only with wing or tail 
loads; the total loads are obtained by adding the fuselage 
loads (which may be assumed to be unaffected by flexibility) 
to the wing or tail loads obtained from the analysis. The 
amount of load carried by a flexible wing and the manner 
of its distribution can consequently be estimated by the 
method presented herein if the contribution of the fuselage 
is known at low dynamic pressure®, that is, for the “rigid- 
wing” case. 

The fuselage has considerable effect on some of the stability 
parameters, although in the case of others, such as C tp , the 
effect is negligible. Other effects that may have to be 
accounted for in calculating stability derivatives are the 
boundary-layer behavior and tip suction. The boundary- 
layer effect may be accounted for by using a section lift- 
curve slope corrected for boundary-layer effects to calculate 
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the angle-of-attack distribution of the flexible wing at the 
flight conditions of interest and then obtaining the lift and 
drag distributions corresponding to that angle-of-attack dis- 
tribution. Lateral tip suction may be important on low- 
aspect-ratio and highly swept wings. Since this suction 
does not affect the lift distribution, it may be taken into 
account by calculating the angle-of-attack distribution of 
the flexible wing and estimating the tip suction corresponding 
to the actual angle of attack at the tip. 

In calculating stability derivatives it is well to keep in 
mind that the method presented in this report is based on. 
a modified strip theory, unless aerodynamic influence- 
coefficient matrices are used . The calculated derivatives may 
therefore, be somewhat in error, particularly if in calculating 
them the moment of a load distribution has to be determined. 
If there is reason to suspect that the modified strip theory 
is inadequate for calculating a given derivative, the deriva- 
tive may be calculated for the rigid-wing case by , a more 
refined method ; the results calculated by the method of this 


report may then be used to correct the accurate rigid-wing 
value for the effect of structural flexibility. 

CONCLUDING REMARKS 

A method has been presented for calculating the aerody- 
namic loading, the divergence speed, and certain stability 
derivatives of swept and unswept wings and tail surfaces of 
arbitrary stiffness. Provisions have been made for using 
either stiffness curves and root-rotation constants or structural 
influence coefficients in the structural part of the analysis. 
Either strip theory with over-all reduction and rounding off 
at the tip or aerodynamic influence coefficients may be used 
for the aerodynamic part of the analysis. Computing forms, 
tables of numerical constants required in the analysis, and 
an illustrative example are included to facilitate calculations 
by means of the method. 

Langley Aeronautical Laboratory, 

National Advisory Committee for Aeronautics, 
Langley Field, Va., December 24, 1948. 


APPENDIX 

SUMMARY OF MATRIX ALGEBRA PERTINENT TO THE ANALYSIS 


For the convenience of the reader unfamiliar with matrix 
method, a summary of matrix definitions and operations is 
presented herein. For a more complete discussion of matrix 
methods the reader is referred to any text on matrices — for 
instance, reference 13. 

DEFINITIONS 

A matrix is a rectangular array of numbers, called ele- 
ments, written in rows and columns. A column matrix 
consists of a single column, a row matrix of a single row. 
A square matrix has as many rows as it has columns. 

The diagonal of a square matrix from the upper left to the 
lower right is called the principal diagonal. A matrix all 
the elements of which are zero except for those on the prin- 
cipal diagonal is called a diagonal matrix. If all of these 
elements are unity, the matrix is termed the unit matrix. 

The transpose of a square matrix is the square matrix 
which results from interchanging the rows and columns in 
the given matrix; it may, consequently, be thought of as 
having been obtained by rotating the given matrix about its 
principal diagonal. 

MATRIX ALGEBRA 

Two matrices can be added or subtracted if both have the 
same number of rows and columns. The addition or sub- 
traction is carried out by adding to or subtracting from each 
element of the first matrix the corresponding element of the 
second matrix. 

A matrix is multiplied by a constant by multiplying each 
element by that constant. 


Two matrices can be multiplied by each other if the 
second has as many rows as the first has columns. Each 
element of the resulting matrix is obtained by multiplying 
the elements in the corresponding row of the first matrix by 
those of the corresponding column of the second matrix in 
the following order: The first element of the row is multi- 
plied by the first element of the column, the second, by the 
second, and so forth. The sum of the products obtained in 
this manner is the value of the element of the product matrix. 
Schematically this process may be illustrated as follows: 


[«] [MJ [»] [.U| 


a b c d e f g 


” . . . . A . . “ 

. ... B . . 

. ... C . . 




X 

. ... D . . 

. ... E . . 

. ... F . . 

. ... G . . 


* 


where 

Q=aA+bB+cC+dD+eE+fF+gG 


It must be emphasized that in multiplying matrices by 
each other their order is of importance. As the two matrices 
under consideration are written, the matrix [m] is said to be 
postmultiplied by the matrix [A/], or the matrix [A/] may be 
said to be premultiplied by the matrix [m]. If the two 
matrices were written in the reverse order and then multi- 
plied according to the foregoing instructions — that is, if the 
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matrix [A/] were postmultiplied by the matrix [m ] — the ele- 
ment of the third row and fifth column of the product 
matrix [21/] [m] would clearly not have the value Q in general; 
nor would, in general, any other element have the value it 
would have if the two matrices were multiplied in the order 
shown. Consequently it is important to observe the order in 
which the matrices are written in the computing instructions. 

MATRIX ITERATION 

The purpose of iterating a square matrix is to determine 
the column matrix or matrices which, if postmultiplied by 
the given square matrix, yield the same column matrix ex- 
cept for a constant multiplier. It is the value or values of 
these multipliers which constitute the desired characteristic 
values of the matrix. 

The iteration is carried out by assuming a “trial'’ column 
(the column shown in table VII (a) is convenient for the pur- 
pose of this analysis) and premultiplying it by the given 
square matrix to yield a “result” column. The elements of 
the result column including the last are divided by the last 
element of the result column and entered as a second trial 
column. The second trial column is then premultiplied by 
the square matrix to yield a second result column. The pro- 
cedure is repeated until the same value (within the desired 
accuracy) is obtained twice in succession for the last element 
of the result matrix. The reciprocal of this value is the de- 
sired (lowest) characteristic value of the matrix, that is, the 
lowest critical value of («z*)n, in the analysis of this report. 

Another way of estimating a first trial column in this 
analysis is to add the elements in each row of the matrix [A], 
enter the six sums in the first result column, and treat them 
as if they had been obtained by multiplying the matrix [A] 
by a first trial column. 

DERIVATION OF THE INTEGRATING MATRICES 

Although familiarity with the derivation of the integrating 
matrices is not essential to the application of the method of 
this report, an outline of the derivation is presented because 
of its general interest. 

The integrating matrices used in this report are based on 
the same concept as Simpson’s rule — replacement of the 
actual function which is to be integrated by parabolic 
segments. If the function y has the values y a - u y X) and 
y n+u respectively, at the equally spaced points * a _ r , *„, and 
the following relations are true for a second-degree 
parabola passed through the three known points: 

y=y*+\ (y.+1-y.-O 

i(^ +1 -2y B +y„_ 1 )(^) 2 (AD 


J X * + ^ydx=(^Ax^y x - 1 +(^Ax^y x +(^Ax^y K+1 (A2) 
j x n * 1 ydx=(— j^Aa:^y»_ l -|-0Aa:^y«-{-^^Aa:^y« + i (A3) 
j x " ^ydx=(j^Ax^y a -i+(^Ax^y x +(—j^Ax^y n+ i (A4) 
J^* + 1 (ar— x x )ydx=(— ^A^iN-i+CfOy.+Qs**^.!.! 

(A5) 

j x * +1 (x —x n ) y dx =(— -^ Sx*)y»-i+^ a* 2 ) y, -f 

(^Ar 2 )?/ B+ I (A6) 

where 

Ax = x tl —x a - l =x a+l —x x (A7) 

The different integrations over the parabolic segments 

may thus be performed by multiplying the given values of 
y by the multiplying factors indicated in equations (A2) 
to (A6). 

Since load distributions at subsonic speeds go to zero with 
infinite slope at the tip and the ordinary second-degree 
parabola - furnishes a poor approximation to such a distri- 
bution, V. M. Falkner has suggested that a curve of the type 

y=A-Mi(l-*) I/8 -M2(l-*) 3/l (AS) 

be passed through the last three points of the load-distribution 
curve at the tip (a?= 1) . On the basis of this approxima- 
tion, relations equivalent to equations (Al) to (A6) may be 
derived. The multiplying factors for the last two segments 
are then based on these equivalent expressions rather than 
those of equations (A2) to (A6). 

The integrating factors of equations (A2) to (A6) may be 
assembled directly into integrating matrices. The matrix 

[/]", for instance, is set up to perform the integration J y dx. 

If at the upper limit *=0.1 and if the ten-point matrix 
(table III (b)) is to be used, the factors 0.04167, 0.06667, 
and —0.00833 may be obtained from equation (A4) since 
*»_i=0, *„=0.1, and A*=0.1; similarly, if for the same case 
the integration is extended to *„ +1 =0.2 as the upper limit, 
the integrating factors 0.03333, 0.13333, and 0.03333 will be 
obtained from equation (A2). These factors constitute the 
second and third rows of the matrix \T \ n since the integra- 
tions are independent of the values of y other than the first 
three, the other values of y are multiplied by zero in these 
two rows. In order to extend the integration to *=0.3 an 
integration is again performed up to *=0.2 and another 
integration, using another parabolic segment, is performed 



924 


REPORT 1000 — NATIONAL ADVISORY COMMITTEE FOR AERONAUTICS 


from a:=0.2 to £=0.3. For the latter integration £„_ 1 =0.2, 
£„=0.3, and A£=0.1, so that equation (A4) again yields the 
factors 0.04167, 0.00667, and —0.00833. The value of y at 
£=0.2 is therefore assigned a multiplying factor of 0.03333 
by the first integration and a factor of 0.04167 by the second, 
or a total factor of 0.07500. The resulting factors are entered 
in the fourth row of the matrix [I]". All other rows are 
obtained in a similar manner. 

The matrices [I] and [/'] are set up to perform the integra- 
tion y dx. The values of the last, row of the ten-point 

matrix [!'] (table I (b)) are obtained from the equivalent of 
equation (A3) for the function given by equation. (A8), 
with £„_i=0.8, £ n =0.9, £,,+!= 1.0, and A£=0.1. Only the 
multiplying factors for the values of y at £=0.8 and £=0.9 
are listed, since the value of y at £=1.0 (the wing tip) is 
assumed to be zero in this analysis, so that its multiplying 
factor is immaterial. The values of the last row but one are 
obtained similarly from the equivalent of equation (A2). 


The values of the row for '^2 = 0.7 are obtained by using 

equation (A3) in the interval £=0.6 to £=0.8 and the equiva- 
lent of equation (A2) in the interval £=0.8 to 1.0. Similarly 


the row for ^ 7 j 2 = ^-^ ^ obtained by combining the results of 

equation (A2) for the interval £=0.6 to 0.8 with the equiva- 
lent. of equation (A2) for the interval £=0.8 to 1.0. All 
other rows are obtained in a similar manner. 

The matrix [7] is obtained in the same manner as the 
matrix [/'], except that equations (A2) and (A3) are used 
at the tip instead of their equivalents. This procedure gives 
rise to a matrix which has one more row and column than the 
matrix [/']. (See tables I and II.) The matrix [Z] performs 
the same operation as the matrix [Z]" except for the direction 
of integration. As a result of this difference the matrix [Z]" 
is essentially a double transpose (first about the principal 
diagonal, then about the other diagonal) of the matrix [I], 
as implied by its symbol. 

The matrices [77] and [II'] are set up to perform the inte- 
gration (£— £o )y dx, where x is the variable of integration 

and £o, the value of x at the lower limit. In applying the 
integrating factors of equations (A2) to CA6) to this integra- 
tion it must be realized that 


/ (*—Xo)y dx= (£„—£») jydx- j-J {x—x n )y dx (A9) 


so that the integrating factors for this integration would be 
obtained by adding (x n —x 0 ) times the factors of equation 
(A2) or (A3) to the factors of equation (A5) or (AG), respec- 
tively, the choice of equations depending on the limits of 
the integration. The factors for the different segments 
(£=0.8 to 1.0, 0.6 to 0.8, and so forth) are then combined for 
any given row (with its given value of £ 0 ) in the manner indi- 
cated for the matrix [/'] to yield the matrix [//']. 

The matrix [C{[ sums up the torques outboard of a given 
point, whereas the matrix [<7 2 ] gives the sum of the moments 
of forces applied outboard of a given point. Neither requires 
any integrations in the sense of equations (A2) to (AC). For 
the six-point method these two matrices are: 


m 


n 

672 

0.2 

04 

0.6 

0.8 

0.9 

1.0 

0 

.1 

1 

1 

1 

1 

1 

0.2 

0 

1 

1 

1 

1 

1 

0.4 

0 

0 

1 

1 

1 

i 

0.6 

0 

0 

0 

1 

1 

1 

OS 

0 

0 

Q 

0 

1 

1 

09 

Q 

a 

0 

0 

0 

1 


[cy 


V 

b '!2 

02 

04 

0.8 

0 .S 

0.9 

1.0 

0 

0.2 

0.4 

OC 

08 

09 

1.0 

0.2 

0 

.2 

.4 

.6 

.7 

.s 

0.4 

0 

0 

.2 

.4 

.6 

.6 

0 6 

0 ■ 

6 

0 

.2 

.3 

.4 

0.8 

0 

0 

0 

0 

.1 

.2 

09 

0 

0 

0 

0 

0 

► ! 


The moment arms which comprise the matrix [C t ] are frac- 
tions of b'/2, so that the matrix must be multiplied by the 
length b'/2 in order to yield actual moments, as staled in 
equation (28). 
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* 

b’t 2 

0 

.2 

.4 

.6 

*8 

.9 

0 

0.06667 

0.26667 

0. 13333 

0.26667 

0.09333 

0.15085 

.2 

-. 01667 

.13333 

.15000 

.26667 

.09333 

.15085 

.4 

0 

0 

.06667 

.26667 

.09333 

. 15035 

.6 

0 

0 

-.01667 

.13333 

.11000 

.15085 

-8 

0 

0 

0 

0 

.02667 

.15085 

.9 

0 

0 

0 

0 

-.01SS6 

.09333 


TABLE I.— VALUES OF THE INTEGRATING MATRICES [/] AND [/'] 

(aj Six-Point Solution 
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TABLE II— VALUES OF THE INTEGRATING MATRICES [//] AND [//'I 



(b) Ten-Point Solution 




(b) Ten-Point Solution 
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TABLE IV.— VALUES OF THE LOAD-CONVERSION MATRIX [C,'] 

(a) Six-Point Solution 


672 

0 

Jl 

.4 


.8 

o 

0.2 

0.08333 

0.13333 

- 0.01667 

0 

0 

0 

A 

-.01667 

.13333 

- OS 333 

0 

0 

0 

A 

0 

0 

.08333 

.13333 

-.01667 

0 

.8 

0 

0 

-.01667 

.13333 

.08333 

0 

.9 

0 

o 

0 

0 

.04553 

.05752 

L 0 

0 

0 

0 

0 

-.01886 

.09333 


(b) Ten-Point Solution 


V 

b '/2 

0 

.1 

.2 

.3 

.4 

.5 

3 

.7 

.8 

.9 

0.1 

0.04166 

0.06667 

' — 0.00833 

0 

0 

0 

0 

0 

0 

0 

.2 

-.00833 

.06667 

.04167 

0 

0 

0 

0 

0 

0 

0 

.3 

0 

0 

.04166 

.06667 

-.00833 

0 

0 

0 

0 

0 

A 

0 

0 

— - OOS 33 

.06667 

.04167 

0 

0 

0 

0 

0 

.5 

0 

0 

0 

0 

.04166 

.06667 

-.00833 

0 

0 

0 

.6 

0 

0 

0 

0 

-.00833 

.06667 

.04167 

0 

0 

0 

.7 

0 

0 

0 

0 

0 

0 

.04166 

.06667 

-.00833 

0 

.8 

Q 

0 

0 

0 

0 

0 

-.00833 

.06667 

.64167 

0 

.9 

0 

0 

0 

0 

0 

0 

0 

0 

.04166 

0 

1.0 

0 

0 

0 

0 

0 

0 

0 

0 

-.00833 

.09333 


TABLE V— VALUES OF THE LOAD-CONVERSION MATRIX [CV] 

(a) Six-Point Solution 


V 

672 

0 

jt 

.4 

.6 

.8 

.9 

0.2 

a 01667 

a 16667 

0.01667 

0 

0 

a 

.4 

-.00833 

.05000 

.11667 

.05000 

-.00833 

0 

.6 

0 

0 

.01667 

.16667 

.01667 

0 

.8 

0 

0 

-.00833 

.05000 

. 0 S 946 

.02035 

.9 

0 

0 

0 

0 

.00631 

.08860 

LO 

Q 

0 

0 

0 

-.01077 

.04190 


(b) Ten-Point Solution 


V 

672 

0 

a 

.2 

.3 

.4 

.5 

.6 

.7 

& 

.9 

0.1 

0.00833 

0.08333 

0.00831 

0 

0 

0 

0 

0 

0 

0 

.2 

-.00417 

.02500 

.05834 

.02500 

— . 00417 

0 

0 

0 

0 

0 

.3 

0 

0 

.00833 

.08333 

.00833 

0 

0 

0 

0 

0 

.4 

0 

0 

-.00417 

.02500 

.05834 

.02500 

-.00417 

0 

0 

0 

.5 

0 

0 

0 

0 

.00833 

.08333 

.00833 

0 

0 

0 

.6 

0 

0 

0 

0 

-.00417 

.02500 

.05834 

.02500 

-.00417 

0 

.7 

0 

0 

0 

0 

0 

0 

.00833 

.08333 

.00835 

0 

.8 

Q 

0 

0 

0 

0 

0 

-.00417 

.02500 

.06629 

.02035 

.9 

0 

0 

0 

0 

0 

0 

0 

0 

.00631 

.08860 

1.0 

0 

0 

0 

« 

0 

0 

0 

0 

-.01077 

.04190 


056646 — SI 60 
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TABLE VI— FORM FOR COMPUTATION OF AEROELASTIC MATRIX — Continued 


(b) Computing Instructions 
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TABLE VI.— FORM FOR COMPUTATION OF AEROELASTIC MATRIX— Continued 

(b) Computing Instructions— Continued 
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TABLE VL— FORM FOR COMPUTATION OF AEROELASTIC MATRIX— Concluded 

(b) Computing Instructions— Concluded 
Supersonic Case 
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TABLE VII.— FORM FOR SOLUTION C 


(a) Divergent 


[A] 

V 

672 

0 

.2 

A 

0 

'■o' 

Q , 

0 

.2 




A 




.6 




.8 _ 




.8 






672 

a) 


op 

0 

0 

0 

0 ; 

.2 

.3000 

1 


.4 

.5000 



.6 

.7000 

BSM 

flgl 

.8 

.9000 



.9 


1.0000 

1.0000 


U)M 

672 

0 ) 

( 2 ) 

( 3 ) 

0 

0 

" 0 

0 - 

.2 




A 




.6 




.8 



■ H 

.9 





(*J *)o ! 
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TABLE VII.— FORM FOR SOLUTION OF AEROELASTIC EQUATION— Continued 


<b) Aerodynamic Loading 
?=, * 9 *= 



Final matrices 
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TABLE IX— COMPUTATION OF AEROELASTIC MATRIX OF EXAMPLE WING (SUBSONIC CASE) 




.04585 

.08216 


mm 

.Q3G68 

.16433 

.03068 

.16433 


.16433 


-.06851 

; 0 

.27401 

0 

.44524 

.40203 



w. 

( <2» r ~ Qa v tan A) 0.1029 


[©] + [©] + [©] 




-.03372 

.23282 

-.05955 

.46567 


.46567 



-.05955 

.46567 

-.05955 

.46567 



mm 

—.00613 

.01777 

.01728 

.01922 

.00503 

.00089 

mm 

-.01173 

.03724 

.04875 

.03089 

.02118 

.02902 


-.01173 

.03724 

.05082 

.17174 

.06273 

.07101 


-.01173 

.03724 

.03938 

.22118 

.10693 

.16293 

.9 'j 

-.01173 

. 03724 

.03938 

.22118 

.10254 

.23190 




-.00199 

.01106 

.02114 

-.00399 

.01543 

.04289 

-.00399 

.01543 

.03755 

-.00399 

.01543 

.03221 

-.00399 

.01543 

.03221 






.42084 I 1.717< 






















































































































































































































































































































TABLE X.— SOLUTION OF AEROELA8TIC 


(a) Divergence 


w 

V 

Fp 

B 

.2 

B 

.6 

.8 

.0 

0 

S 

0 

wm 

0 

0 

0 

.2 

-.0041 

.0007 

0072 

-.0028 

— , 0800 

-.0600 

n| 

-.0077 

.0218 

,0060 

1259 

-.0787 

-, 1541 

ini 


.0218 

,0188 

1122 

-.1027 

-.2008 

, 8 . 

-.0077 

.0218 

.0072 

-. 08)3 

-.0802 

-.3204 

.0 

-.0077 

.0218 

.0072 

-.0818 

-.0848 

-.2902 


M 

V 

672 

(» 

( 2 ) 

( 3 ) 

( 4 ) 

( 8 ), 

( 0 ) 

0 

0 

0 

0 

0 

0 

0 

.2 

.3000 

.3116 

.8449 

,3480 



flEKj 

.6000 

.7811 

.7840 

.7870 



S3 

.7000 

1.0280 

1. 0626 

IjUfljJSl 



.8 

.9000 

1.0776 

1.0714 

1.0718 



.9 

1.0000 

1, 0000 

1.0000 

1.0000 

1.0000 

1.0000 



0 

,2 

,4 

"To 

.8 

.8 


MM 

*1 

T /2 

(1) 

(2) 

(3) 

(4) 

(6) 

(0) 

0 

wm 

0 

0 

0 

0 

0 

,2 

-.1286 

-.1501 

-. 1670 




,4 

l 

CO 

o 

CO 

-.8551 

-.3607 




.6 

-.4240 

-.4704 

-.4770 




.8 

-.4448 

-.4849 

-.4862 




.9 

-.4128 

-.4620 

-.4529 

-.4529 




- 2.208 



EQUATION FOR EXAMPLE WING (SUBSONIC CASE) 

(b) Aerodynamic Loading 
-®-— 0.25 ?- K 5 *- 0 ,. 1.12 


[U 1 -wM] 



0.0007 

0.2418 

0. 1081 

0.1909 

0.0881 


L^'l] Ml 


.0433 

0. 1140 

0. 0466 

0.0792 


1®J{S 1-0,8524 

l®ll“i" 

L®J 1 «)“ 

-0.2730 

1®J(«| - 0,1681 

1 ®IH- 






C/,„-0.740C z , a 

o B M„-o.mc La 

<7^-0,0855(7^ 

5^2"°' 41S 


eo 

co 


















































































































































































































































